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Abstract 

New checkable criteria for persistence of chemical reaction networks are proposed, which 
extend and complement those obtained by the authors in previous work. The new results 
allow the consideration of reaction rates which are time-varying, thus incorporating the 
effects of external signals, and also relax the assumption of existence of global conservation 
laws, thus allowing for inflows (production) and outflows (degradation). For time-invariant 
networks parameter-dependent conditions for persistence of certain classes of networks are 
provided. As an illustration, two networks arising in the systems biology literature are 
analyzed, namely a hypoxia and an apoptosis network. 



1 Introduction 

For differential equations evolving in Euclidean space, "persistence" is the property that all 
solutions starting in the positive orthant do not approach the boundary of the orthant. In- 
terpreted for chemical reactions and population models, this translates into a "non-extinction 
property" that states that no species will tend to be completely eliminated in the course of the 
reaction, provided that every species was present at the start of the reaction. 

In the previous work [2] , we presented criteria for checking persistence in closed chemical reac- 
tion networks, couched in the language of graph theory and Petri nets. One of the main results 
was that a time-invariant, conservative chemical network is persistent provided that each siphon 
contains the support of a P-seminow, regardless of the reaction kinetics underlying the chemical 
reactions, or the values of parameters such as rate constants. 

In the present paper, we extend the previous results in several directions: (1) Kinetic coefficients 
are now allowed to be time-varying. (2) No conservation assumption is made (so solutions are 
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potentially unbounded). (3) The case in which there are critical siphons (i.e., do not contain 
the support of any P-semiflow) is studied, and a sufficient condition for persistence is provided 
in that case. This latter condition is parameter-dependent. 

The motivation for considering time-dependent coefficients is that these may be used to repre- 
sent the effect of external inputs to the network, while not making a conservation hypothesis 
allows the consideration of inflows and outflows, or production and degradation processes. Fi- 
nally, critical siphons arise in many examples. 

We provide two examples of applications of the new results. The first one analyzes a model of the 
common core subsystem responsible for the hypoxia control network in C.elegans, Drosophila, 
and humans. Hypoxia (deprivation of adequate oxygen supply) results in the expression of 
specific genes in response to stress caused by low concentration of available oxygen. This 
particular example was picked for two reasons. First of all, viewing oxygen concentration as 
an external input gives rise to a network with time-dependent kinetic coefficients. Second, 
in this model there are no conservation laws that guarantee boundedness of solutions, and so 
the same example serves to illustrate the role of the new concepts of conditional persistence 
introduced in this work. The second example is an apoptosis (programmed cell death) network. 
We characterize persistence using our result for critical siphons. 



2 Background on chemical reaction networks 

A chemical reaction network ("CRN", for short) is a list of chemical reactions Ri, taking 
place among species Sj, where the indices i and j take values in 7Z := {1,2, ...,n r } and 
S := {1, 2, . . . n s } respectively. Individual reactions are then denoted as follows: 

Ri : J^aijSj ^J^PijSj (1) 
jes jes 

where the ccjj and are nonnegative integers called the stoichiometry coefficients. The species 
j on the left-hand side for which > are called reactants and the ones on the right-hand side 
for which (3 > the products, of the reaction. Informally speaking, the forward arrow means 
that the transformation of reactants into products only happens in the direction of the arrow. 
If also the converse transformation occurs, then, the reaction is reversible and we need to also 
list its inverse in the network as a separate reaction. Sometimes, for convenience of notation, we 
will associate to a reaction the two integer column vectors ccj and fa, whose entries are defined 
by the stoichiometry coefficients. It is worth pointing out that we allow chemical reactions in 
which both the right and left hand sides are actually empty (though not at the same time). 
This case corresponds, from a physical point of view, to inflows and outflows of the chemical 
reaction. 

As usually done, we arrange the stoichiometry coefficients into an n s x n r matrix, called the 
stoichiometry matrix T, defined as follows: 

F\ j i = Pi j -a ij , (2) 

for all i € 1Z and all j G <S (notice the reversal of indices). This will be later used in order 
to synthetically write the differential equation associated to a given chemical network. Notice 
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that we allow T to have columns which differ only by their sign; this happens when there are 
reversible reactions in the network. 

We discuss, next, how the speed of reactions is affected by the concentrations of the different 
species. Each chemical reaction takes place continuously in time, at its own rate, which is 
assumed to be only a function of the concentration of the species taking part in it. In order to 
make this more precise, we define the vector S = [Si, S2, ■ ■ ■ S na ]' of species concentrations and, 
as a function of it, the vector of reaction rates 

R(S,t) := [R 1 (S,t),R 2 (S,t),...R nr {S,t)Y 

where t G [0, 00) denotes time. Notice that we explicitly allow time-dependence, as we wish 
to to consider the effect of external inputs to the system. Such inputs may represent chemical 
species which are not explicitly considered as part of the state variables but which, nevertheless, 
influence the reaction rates. 

Some mild uniformity requirements are needed for technical reasons as far as time time- 
dependence is concerned. We assume, in particular, that for alH G 1Z for all S and alH > 

Ri(0, t)=0 and R^S) < Ri(S, t) < R^S), (3) 

where the R i (S),Ri(S) are non-negative, continuous functions of S, satisfying the following 
monotonicity constraint: 

S^ Ri S ^Ri{S) >Ri{S) (4) 

for alH G 1Z (and a similarly for R^), where the notation S 2>^ 4 S means that we have a strict 
inequality Sj > Sj whenever species j is a reactant in reaction i. (We also write, more generally, 
S S> S for any two vectors of species concentrations, if Sj > Sj for all j = 1, . . . , n s .) 

Furthermore, we assume standard regularity assumptions of Ri(S,t) in order to ensure local 
existence and uniqueness of solutions. 

A special form of reaction rates are mass-action kinetics, which correspond to the following 
expression: 

n s 

Ri(S, t) = ki(t) Yl Sj' 3 for alH = 1, . . . ,n r 

3=1 

(interpreting S* = 1 for all S), that is, the speed of each reaction is proportional to the 
concentration of its reagents. Notice that we allow a time- varying kinetic rate ki(t), which 
may account for the effect of external species not explicitly included in the network under 
consideration. In the case of mass-action kinetics, a uniform lower and upper bound on Ri(S, t) 
exists if and only if there exist constants k\ ni > and fc* up > such that 

4rf<*(t)<C P vt>o. 

With the above notations, the chemical reaction network can be described by the following 
system of differential equations: 

S(t) =TR(S(t),t). (5) 
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where S = S(t) evolves in M™ and represents the vector of all species concentrations at time 
t, and r is the stoichiometry matrix. For systems with mass-action kinetics the following 
alternative expression is valid: 

S(t) = J2(Pi-aiMt)S ai (t) (6) 

ieTZ 

where f3i = is the column vector co\((3n, • • • , A'n s )> a i = a i* is the column vector col(aji, . . . , aj n J, 
and 5 7 = Sj 1 . . . SZ" S for any nonnegative vector 7 = (71 , . . . , 7 ns ) . 

It is straightforward to verify that the positive orthant is positively invariant for system ([5]). 
Moreover, for each So E M" the affine subspace defined by: 

So + im[r] 

is also invariant, regardless of the specific expression of reaction rates; its intersection with the 
positive orthant (which is therefore a forward invariant set) is called the stoichiometry class of 
S . 

3 Petri nets and structural invariants 

In stating our results, we will employ some terminology borrowed from the graph theory, and 
specifically Petri nets. Although arising from the study of discrete processes, Petri nets provide 
a useful language and graphical representation for CRN's, and a number of structural and 
analytical tools developed for them can be easily adapted to the continuous context of chemical 
reactions. In what follows, we associate to a CRN a bipartite directed graph (i.e., a directed 
graph with two types of nodes) with weighted edges, called the species-reaction Petri net, or 
SR-net for short. Mathematically, this is a quadruple 

(V S ,V R ,E,W), 

where Vs is a finite set of nodes, each one associated to a species, Vr similarly is a finite set of 
nodes (disjoint from Vs) corresponding to reactions, and E is a set of edges as described below. 
(We often write 5 or Vs interchangeably, or R instead of Vr, by identifying species or reactions 
with their respective indexes; the context should make the meaning clear.) The set of all nodes 
is also denoted by V = Vr U Vs- 

The edge set E C V x V is defined as follows. Whenever a certain reaction Ri belongs to the 
CRN: 

^aijSj -► ^PijSj, (7) 

we draw an edge from Sj E Vs to Ri E Vr for all Sj's such that cty > 0. That is, (Sj, Ri) E E 
iff dij > 0, and we say in this case that Ri is an output reaction for Sj. Similarly, we draw an 
edge from Ri E Vr to every Sj E Vs such that /3y > 0. That is, (Ri, Sj) E E whenever /3jj > 0, 
and we say in this case that Ri is an input reaction for Sj. 

Notice that edges only connect species to reactions and vice versa, but never connect two species 
or two reactions. 
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More generally, given a nonempty subset E C S of species, we say that a reaction Pj is an 
output (input) reaction for E if it is an output (input) reaction to some species of E. 

The last element to fully define the Petri net is the function W : E — > N, which associates to 
each edge a positive integer according to the rule: 

W(S j} Ri) = onj and W(R U S 3 ) = fly . 

The stoichiometry matrix T, previously introduced, is usually referred to as incidence matrix 
in the Petri Net literature. 

Several other definitions which are commonly used in that context will be of interest in the 
following. We say that a row or column vector v is non-negative, and we denote it by f ^ 
if it is so entry- wise. We write v y if v y and v ^ 0. A stronger notion is instead d>0, 
which indicates Vi > for all i. 

Definition 3.1 A P-semiflow is a row vector c y such that cT = 0. The support of a 
P-semiflow is the set of indexes {i € Vs : Cj > 0}. 

Using the fact that the entries of T are integers, it is easy to show that, given any P-semiflow 
c, there is always a P-semiflow with integer components which has the same support as c. 

Definition 3.2 A nonempty subset E C S of species is stoichiometrically constrained if there is 
a P-semiflow whose support is included in E. When E = 5, that is, if there is some P-semiflow 
c S> 0, we simply say that the CRN (or the corresponding Petri net) is stoichiometrically 
constrained. 

P-semiflows for the system ([5]) correspond to non-negative linear first integrals, that is, linear 
functions S i— > cS such that (d/dt)cS(t) = along all solutions of ©. In particular, a Petri 
net is stoichiometrically constrained if and only if there is a positive linear conserved quantity 
for the system. (Petri net theory views Petri nets as "token-passing" systems, and, in that 
context, P-semiflows, also called place-invariants, amount to conservation relations for the 
"place markings" of the network, that show how many tokens there are in each "place," the 
nodes associated to species in SR-nets. We do not make use of this interpretation in this paper.) 

Definition 3.3 A T-semiflow is a column vector v y such that Tv = 0. 

Once again, one can assume without loss of generality that such a v has integer entries. 

Definition 3.4 A nonempty subset A C TZ of reactions is consistent if there is a T-semiflow 
whose support includes A. When A = TZ, we also say that the CRN, or its associated Petri net, 
is consistent. 

The notion of T-semiflow corresponds to the existence of a collection of positive reaction rates 
which do not produce any variation in the concentrations of the species. In other words, v can 
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be viewed as a set of fluxes that is in equilibrium [9j. (In Petri net theory, the terminology is 
"T-invariant," and the fluxes are flows of tokens.) 

A vector v = (0, 0, . . . , 0, 1, 0, . . . , 0) with a "I" in the ith position and 0's elsewhere represents 
the ith reaction; thus we may label such a unit vector as "iV . With this notational convention, 
the following fact holds. Suppose that and Ri are reactions that are reverses of each other, 
that is, ctkj = Pej and fl^j = oigj for every species j G S. Then, the vector R^ + Rg is a 
T-semiflow, becaus the kth and ith columns of T are opposites of each other. In chemical 
network models of biological systems, it is common for several of the reactions to be considered 
as reversible. This gives rise to many such "trivial" T-semiflows. 

Definition 3.5 A nonempty set £ C Vs is called a siphon if each input reaction for £ is also an 
output reaction for £. A siphon is minimal if it does not contain (strictly) any other siphons. 

For later use, we associate a particular set to a siphon £ as follows: 

L s = {x G M> s \xi = i G £}. 

The set is therefore characterized as the set of concentration vectors whose entries are zero 
if (and only if) the corresponding chemical species are in the siphon £. 

4 Persistence and consistency 

Our main interest is the study of persistence: when do species remain nonzero, if they start 
nonzero? We will study two variants of this concept, and will provide a necessary characteri- 
zation for one and a sufficient characterization for the other. 

Definition 4.1 A nonempty subset £ C S of species is conditionally persistent (respectively, 
bounded-persistent) if there exists a bounded solution S{ ) with S(0) 3> such that 

liminf S,-(t) > VjG£ (8) 

t — >oo 

(respectively, if this property holds for all bounded solutions). When £ = <S, we say simply 
that the corresponding CRN is conditionally persistent or bounded-persistent respectively. 

Note that in case £ = S condition ([8]) amounts to the requirement that the omega-limit set 
uj(S(0)) should not intersect the boundary of the main orthant. 

The following result generalizes Theorem 1 in j!2j to systems with time- varying rates, and 
provides a necessary condition for persistence. It is proved in Section G3 

Theorem 1 Every conditionally persistent CRN is consistent. 

The following result gives a sufficient condition for persistence, and generalizes Theorem 2 in [2] 
to systems with time- varying rates: It is proved in Section [5j 

Theorem 2 If a CNR has the property that every siphon is stoichiometrically constrained, then 
it is bounded-persistent. 
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5 Key technical results 



Definition 5.1 A nonempty subset A C TZ of reactions is conditionally persistent if there exists 
a bounded solution 5(-) with 5(0) ^> such that 

\unm.fRi(S(t),t) > Vi € A . 

The first key technical fact that we need is as follows; it is proved in Section [3 

Theorem 3 Every conditionally persistent subset A of reactions is consistent. 

Definition 5.2 A nonempty subset E C S of species is an extinction set if there exists a 
bounded solution 5(-) with 5(0) 3> and a sequence t n — > oo such that 

lim 5,(i„) = O jGS. (9) 

n— >oo 

Equivalently, E is an extinction set if and only if Ls f"|u;(5(0)) 7^ for some bounded solution 

The second key technical fact, proved in Section [6J is as follows. 

Theorem 4 Every extinction set of species is a siphon. 

Lemma 5.3 An extinction set cannot be stoichiometrically constrained. 

Proof. Let E be an extinction set, and pick a bounded solution 5(-) as in the definition of 
extinction. Suppose that E C S is stoichiometrically constrained. Let c be a P-semiflow whose 
support is included in E. Since c is a P-semiflow, cS(t) = c5(0) > for all t > (the last 
inequality because 5(0) S> and c > 0). Since the support of c is a subset of E, it follows that 
c5(i) = CjSj(t), with the sum only over the indices j G E. Thus ([9]) cannot hold. I 

We defined what it means for a reaction Ri to be an output reaction for a species Sj, namely 
that Sj should be a reactant of Ri. More generally, we use the following concept. 

Definition 5.4 Consider a nonempty subset E C S of species. A reaction Ri is said to be a 
sink for S if the set of reactants of Ri is a subset of E. The set of all sinks for E is denoted as 
A(E). 

Lemma 5.5 If E is conditionally persistent, then A(E) is conditionally persistent. 

Proof. Suppose that we have a bounded solution such that liminfi_^oo Sj(t) > s > for every 
j € E. Since the solution S(-) is bounded, this means its closure is a compact subset K of 
the (closed) positive orthant. Pick any sink Ri for E. By Property Ri(S) > for all 
5 6 K. Therefore liminf^ 0O i?j(5(t),t) > for this same trajectory, which proves that A(E) 
is conditionally persistent. I 
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Corollary 5.6 If S is conditionally persistent, then A(S) is consistent. 

Proof. This follows immediately from Lemma 15.51 and Theorem [3l I 

This completes the proof of Theorem[TJ because the hypothesis of the Theorem says that £ = S 
is conditionally persistent. By Corollary 15.61 A(«S) = 7Z is consistent, which means that the 
CRN is consistent, as claimed. 

To prove Theorem [2j we observe: 

Lemma 5.7 If a CRN is not bounded-persistent, then there is some extinction set. 

Proof. Suppose that there is some bounded solution S(-) with S(0) S> 0, some species jo € S, 
and some sequence t n — ► oo, such that linv^oo Sj(t n ) = 0. For this solution, and for this same 
sequence {t n }, let S be defined as the set of species j G S such that linv^oo Sj(t n ) = 0. Since 
jo £ X, £ is nonempty, and it is an extinction set by definition. I 

Now Theorem [2] follows from: 

Corollary 5.8 If a CRN is not bounded-persistent, then there is a non-stoichiometrically con- 
strained siphon. 

Proof. Assume that the given CNR is not bounded-persistent. By Lemma 15.71 there is an 
extinction set S. By Theorem HI £ is a siphon. By Lemma 15. 3| £ is not stoichiometrically 
constrained. I 



6 Proof of Theorem [4] 

Let the nonempty subset S C S of species be an extinction set. Pick a bounded solution S(-) 
with S(0) » such that Ls P|o;(5(0)) ^ . We need to prove that E is a siphon. 

Assume that y G Ls nw(5(0)) but that S is not a siphon. Hence, there exists a species Sj G £ 
so that for at least one of its input reactions Rk and all of RkS reactant species Si, it holds 
yi > 0. By Property j3|), we have that 

Rk(y,t) >R k (y) = r>o 

for some positive value f and all £ > to- 

Therefore, since all output reactions of Sj have zero rate at y (no matter what the value of t 
is), and at least some incoming reaction is strictly positive, it follows by continuity of each of 
the i?fc(5*)'s that there is some e > 0, so that: 

Sj(t) = [TR{z{t),t)] j >f/2 
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whenever z(t) E B £ (y) = {z y : |z — y\ < e} and i > to- Now, using the uniform upper 
bound R(S) = (Ri(S), . . . , R r (S)) and its continuity, we know that there exists M > so that 
\TR(z,t)\ < M for all z as before. Hence, 



\S(t b ) - S(t a )\ 



TR(S(t),t)dt 



< (t b -t a )M (10) 



whenever S(t) € B £ {y) for i G [ta,fy>]- 

Assume without loss of generality (choosing a smaller e if necessary) that e is such that S(0) 
B £ (y). Consider now any partial trajectory crossing the boundary of B e (y) at time t £ , and 
hitting the boundary of B £ /2{y) at time t £ /2, where t s /2 is picked as the first time after t £ when 
this happens. Notice that such a partial trajectory exists, because S(0) B £ (y) and by our 
assumption that y G uj(S(0)). 

Since Sj(t) > f/2 for all t > to whenever z(£) belongs to B e (y), it follows that necessarily we 
must exit B £ (y) an infinite number of times, hence infinitely many such partial trajectories 
exist. 

By the estimate in (fTO]) . the time it takes to get from the boundary of B £ (y) to B £ /2(y) is at 
least e/2M. Moreover, since Sj(t) > f/2 we have: 



Sj(t £ ) + er/4M > ef/AM. 



Sj(t £/2 ) = s^ + ^fs^dt > 



Obviously, for t > t e / 2 , and as long as S(t) € B £ (y), we also have: Sj(t) > Sj(t £ / 2 ) > er/AM. 
This shows indeed y ^ cj(5(0)), contradicting our hypothesis. Hence, £ must be a siphon. I 



7 Proof of Theorem [3] 

Suppose that A C TZ is a conditionally persistent set of reactions, and pick a bounded solution 
S(-) with S*(0) S> such that liminft-,.^ Ri(S(t), t) > for each j G A. We need to show that 
A is consistent. 

Pick an arbitrary t > to. Clearly: 

S(t,t ,S ) - S = f S{t)dt = T I R{S{t),t)dt. (11) 

J to J t(, 

Since S(t,to,So) is bounded, so are R(S(t),t) and its average 



t 

R(S(t),t)dt. 

to 



Hence, there exists a sequence t n — > +oo such that 



1 



R(S(t),t)dt 
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also admits a limit R y as n — ► +00. Now taking limits along this subsequence in both sides 
of (fTT|) . after dividing by t yields: 

0= lim ^(S(t n ,t ,S ) -S ) = TR, (12) 

n— >+oo t n 



Moreover, for all i G A we have 



liminfi? i (S'(t),t) = r< > 

t— >+oo 



and hence there exists T > so that for all t >T, Ri(S(t),t) > rj/2. Letting r be the minimum 
of the fj's, we conclude that: 

5i = lim n ^ +00 i jj n flj- (5(t) , t) dt 
> lim^(C-to)§ = § > 0. 

So, R is a T-semiflow t> whose support contains A. I 



8 An example: hypoxia network 

As discussed in the Introduction, we analyze a model of the hypoxia control network. Starting 
from the model given in [7J for the core subsystem of the hypoxia control network in C. elegans, 
Drosophila, and humans, with 23 species and 32 reactions, the authors of [8] picked a subsys- 
tem consisting of 13 species and 19 reactions which constitute the key components explaining 
experimentally observed behaviors. We analyze this simplified model. 

One of the species, Si, which represents the transcription factor HIFa, is subject to production 
and degradation (or, in formal terms, "inflows" and "outflows"). 

The reactions are as follows: 






- Si 





Si + S 7 


<-> S2 




S2 + S4 






Si + S 6 


*-> s 8 


— > Sq + Sn 


S4 + Sio 


*-> Sg 




5*10 


<-> Sn + 57 




S2 + Se 


<-> S 5 


— > 56 + 5io 


Sn + S12 


*-> S13 


—* S12, 



(13) 



where the meaning of the various biochemical species is in Table [TJ 

External oxygen affects the dynamics of the system by scaling the rate constants for the reactions 
Sg — * Sq + Sn and S$ — > 56 + 5io- Mathematically, this means that ki(t) is proportional to the 
oxygen concentration (and hence is potentially time-dependent) for each of these two reactions. 

According to [7], when the oxygen level falls below a critical value, a sharp rise in HIFq is 
observed, while this protein is undetectable if the oxygen level is above the critical value. The 
modeling effort in [7J and the analysis of the model in [8] were aimed at understanding this 
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51 


HIFa 


S2 


HIFa:ARNT 


S3 


HIFa:ARNT:HRE 


54 


HRE 


55 


HIFa:ARNT:PHD 


56 


PHD 


57 


ARNT 


58 


TTTTl TAT T I "\ 

HIFa:PHD 


59 


HIFa OH:ARNT:HRE 


510 


HIFa OH:ARNT 


511 


HIFa OH 


512 


VHL 


513 


HIFa OH:VHL 



Table 1: The various species in the hypoxia network f)13[) . 



switch-like behavior. For simplicity, the oxygen level was kept constant, and the model was 
investigated over a range of values of this constant. Here we will investigate the persistence 
properties of this network under the assumption that oxygen levels are time-varying. 

The associated reaction network, represented as a Petri Net, is shown in Fig. [TJ This network 
admits 4 linearly independent P-semiflows, which are associated to the following conservation 
laws: 

5i2 + 5i3 = consti 
5 3 + 5 4 + Sg = const 2 , .x 

55 + 5 6 + 5 8 = const 3 
52 + 5 3 + 5 5 + 5 7 + Sg + 5io = const4 

and it is clearly not conservative, due to the presence of outflows and inflows (equation for 
Si). Not only are there no strictly positive conservation laws vT = 0, but there are not even 
"decreasing" semiflows satisfying vT < which could be used as Lyapunov functions in order 
to establish boundedness of solutions. In addition, some kinetic rates are allowed to be time- 
varying. Thus, the techniques from [2] cannot be applied to study persistence. 

We consider next the possible T-semiflows. There are several "trivial" ones, corresponding to 
the reversible reactions: i? 3 + R4, R§ + Rq, Rj + R$, Rio + R\\, R12 + Ri3,Ru + R15 and 
Ru + Rig. 

In addition to these, one can find 3 non-trivial independent T-semiflows: 

i?i + R2 

Ri + R 7 + R 9 + R u + Ris (15) 
R\ + R3 + R15 + Rie + R12 + Rn + Ris ■ 

Since every reaction appears in at least some T-semiflow, the sum of the semiflows shown is 
also a semiflow which is strictly positive, and we can conclude that the network is consistent. 

Thus, the necessary condition for persistence in Theorem [3] is satisfied. This does not quite 
prove conditional persistence, but shows that the property is not ruled out by the structure of 
the network. 
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Figure 1: Hypoxia network 



Next, we find a set of minimal siphons: 

{5*12, S13} 

{S3, S4, 69} , . 

{S2, S3, S5, 1S7, Sg, Sio}- 

Notice that all of them coincide with the support of some P-semiflow, hence every siphon 
is stoichiometrically constrained. We conclude that the network is bounded-persistent, by 
Theorem[2j The next step is the investigation of which variables have the potential for becoming 
unbounded. An algorithm developed for this purpose is illustrated in [1J. It is based on a linear 
time-varying embedding of individual species equations and it carries out a consistency check 
in order to verify which scenarios are compatible with the topology of the network, assuming 
(for instance) mass-action kinetics. Scenarios are described by labeling each species with a 
symbol in {0, l,u>} depending on its asymptotic behaviour, namely converging to 0, bounded 
and bounded away of or diverging to infinity. While such a classification of behaviours does 
not cover all the potential asymptotic dynamics of general systems, it appears to be, in practical 
situations, a fairly mild restriction. 

Running the algorithm on the hypoxia network yields 3 potential scenarios for the asymptotic 
behaviour of the network: 
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Figure 2: Reduced Petri Net and P-decreasing vector support 



Scenario 
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Moreover, see pQ further details, asymptotic invariant vector analysis allows to discard Scenario 
III, see Fig. [2]). Indeed, the linear function £2 + 5*3 + £5 is associated to a P-decreasing 
vector of the reduced net obtained by removing reaction i?3, which, according to the labeling is 
asymptotically switched off. Hence, in this scenario, the quantity S2+S3+S5 gets asymptotically 
dissipated with a strictly positive rate; a clear contradiction. 

Indeed, simulations showed that both scenarios I and II are possible for different values of the 
kinetic constants. Notice that persistence is violated in scenario II as several species (namely 
5*2, S3, £5, £7, S12 ) vanish asymptotically. This is not in contrast with our theoretical develop- 
ments as species Su in this case gets accumulated and diverges to infinity, thus violating the 
boundedness assumption which is crucial to Theorem [2j This also shows that extinction sets 
for unbounded solutions need not be siphons. 
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9 Cascade decompositions and siphons 



One of the main results in [2] states that a time-invariant, conservative CRN is persistent if each 
siphon contains the support of a P-semiflow, regardless of the reaction kinetics underlying the 
chemical reactions, or the values of parameters such as rate constants. In the previous sections 
we have extended this result to CRN's with time-varying reaction kinetics which may have 
inflows and outflows. In this section, we return to the time-independent case, and we explore 
the scenario of a CRN that has a siphon which does not contain the support of a P-semiflow 
(using the terminology of [2], we say that this siphon is critical). We will see that a simple 
sufficient condition based on a linearization argument can be formulated that still guarantees 
persistence of the CRN. The linearization condition is related to, though different from, other 
conditions imposed on chemical reaction networks in order to obtain persistence, notably that 
given in Section 7.2.1 in [4j. 

Let E be a critical siphon. We partition the state S according to species which do not belong 
to E, denoted by (, and species which belong to E, denoted by a. As shown in [2j, siphons 
have the forward invariance property stated below: 

S eL E ^ S(t,S Q ) € L E Vt > 
If the system equations of the CRN are 









<7 







then the forward invariance property implies the following condition: 

/ CT (C,0)=0, VC^O; 



(17) 



(18) 



Consequently, the dynamics on the closure of Ly, (which is a forward invariant orthant of lower 
dimension, as proved in [2]) is completely determined by the equations: 



C(i) = / c (C(i),o). 



(19) 



If co(So) n L-£ is non-empty for some So 3> 0, then it is also a forward invariant set, and so is its 
closure, and the dynamics on the latter set are described by the evolution of the £ components, 
as given by (I19p . Checking non-emptiness of such a set may in general be a rather challenging 
task. However, in some cases it turns out that this is indeed possible by analysis of a reduced 
variational equation. Since £ in (118p is arbitrary, we have the following further result: 



o, 



vc t o. 



(20) 



cr=0 



Therefore, the variational equation associated to solutions corresponding to any initial condition 
S(0) = [C(0)',0]', takes the following block triangular structure: 







6a(t) _ 





^(C(*),o) 3£(C(f),o) 
o &(C(t),o) 



sm 

5a {t) 



(21) 
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These observations suggest that persistence may be understood by examining the stability 
properties of a reduced order variational equation: 



5a(t) = ^(C(t),0)5a(t), (22) 

whose dynamics may determine whether the part of the boundary where a = is repelling or 
attracting to interior solutions. We claim that (|22p is a time-varying positive system because 
the matrix 

|^(C,0) (23) 

is Metzler for all C >z 0. This will follow from the forward invariance of the non-negative orthant 
for system (I17p . Indeed, a standard first-order Taylor expansion yields: 



MM = fa((,0) + ^r(Ma + o(\a 
= ^( C ,0)a + O (H). 



a/„,^_r.„_^ ( 24 ) 



If [-^r(C> would be negative for some integers i ^ j and some £ and letting a = eej, (e.j is 
the j-th element of the canonical basis of the Euclidean space of compatible dimension), yields, 
thanks to (|24j) . /,j(^,ee,)j < for all sufficiently small e > 0. This violates forward invariance 
of the positive orthant (since, clearly, [eej]i = 0). 

Assume that closed invariant sets of (|19p are equilibria, let us denote them by £ e . In view of 
what has been said so far, one is lead naturally to consider the following implications concerning 
persistence : 

1. ^r(CejO) is Hurwitz =4> CRN is not persistent. 

2. 7^(Ce>0) has a positive dominant eigenvalue and is irreducible => Q e £ u(So), 

This suggests that persistence can sometimes be analyzed by simple linearization techniques 
around boundary equilibria. Item 1. of the above claim is obvious, as the existence of a boundary 
equilibrium point with a non-trivial stable manifold clearly violates persistence of a CRN. 

We proceed next to a formal statement and proof of item 2. The proof relies on a few auxiliary 
results which are deferred to the Appendix. 

Theorem 5 Consider a time-invariant CRN whose associated Petri-Net is conservative. Let 
£ be a minimal siphon and let the state S(t) be partitioned accordingly: [C(t), o~(t)]. Let £ e be an 
equilibrium of $19\). globally asymptotically stable relative to its stoichiometry class. Assume, 
moreover, that 

9 fa- ( 9 fa \ 

— — (CejO) be irreducible and Xpp ( -g — (( e ,0)j > (25) 

Then, for any initial condition So = [Co, o"o] in the interior of the positive orthant and denoting 
by lo(So) the uj -limit set of the corresponding solution S(t,So), there holds that S e = [Ce>0] ^ 

"(So). ■ 
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Proof. The proof is by contradiction. Let So in the interior of the positive orthant be such that 
S(t n ,So) — > S e along some increasing sequence t n — > +00 as n — > +00. Let c > be a left 
eigenvector of ^-(Ce,0): 

c ^ (Ce ,0) = A„(^( Cs ,0))c 
and define the following function: 

V(S) := [0,c]S. 

Then V(S e ) = (since S e is of the form [*,0]'), and V(S) > for all S G Q n int(R| ), where 
Q is a neighborhood of S e chosen as in the proof of Lemma IA.ll applied to the vector function 
/<t(C>°")- Furthermore, we assume without loss of generality that Q is relatively open in R> . 
Returning to the solution starting at So, we claim the following: 

There is some p G uj(Sq) n Le wif/i p 7^ S e . 

By passing to a subsequence if necessary we assume that S(t n , So) G Q for all n. There is some 
sufficiently large n* such that for all n > n* , we can define 

r n := sup{t 1 1 < t n and So) ^ Q}, 

which can be thought of as the most recent entry time (into Q) before time t n . Notice that n* 
exists, so that the definition of the T n 's makes sense. Indeed, if n* did not exist, then S(t,So) 
would belong to the open set Qnint(R> ) for all t > 0. In that set, the function V is increasing 
along solutions as we have remarked earlier, and thus limt_* +00 V(S(t, So)) > V(So) > 0, 
contradicting that liminft^ +00 V(S(t, So)) = since S(t n ,So) — > S e and V(S e ) = 0. A similar 
argument shows that t„ — * +00 as n — > +00. 

For all n > n* and t G (r n , i n ), we have that 5(i, 5o) G Q fl int(R> ) and hence V is increasing 
along this part of the solution. In particular, for all n > n*: 

0<V(S(T n ,S )) <V(S(t n ,S )), 

where the quantity on the right tends to zero as n — > +00. Consequently, 

V(S(r n , So)) — > as n — > +00. 

By the definition of V, it follows that V(S(r n , So)) = cp n , where S(r n ,So) = [*,p n ]' for all 
n > n* . Since c ^> 0, it follows that p n — > as n — > +00. 

By continuity of solutions, the sequence S(r n , So) belongs to Q c , the complement of Q, which is 
a closed set. Since the sequence S(r n , So) is bounded, we can pass to a subsequence if necessary, 
and assume that S(r n , So) — > p G Q c as n — > +00. By the previous paragraph, p must belong 
to Ls, and p clearly belongs to u>(So) as well. Finally, since p G Q c , there holds that p ^ S e , 
which establishes the claim. 

On the other hand, let C := lo(So) ncl(Ls). By definition C is a compact set. We wish to show 
that C is invariant. First, any solution starting in C exists for all t G R since it belongs to the 
compact invariant set uj{Sq). Moreover, by (fT8|) and denoting by a(t) the u-component of any 
such solution at time t, it holds a(t) = for all t G R, thus showing that such solutions remain 
in cl(-Ls) for all t G R. It follows that C is indeed invariant. Hence, Corollary IA.31 applied with 
C equal to the closure of co(So) n Ls, and xo = S e , implies that C = {S e }, which contradicts 
the fact that p belongs C. I 
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Figure 3: Apoptosis Network 



10 Example: Apoptosis Regulation Pathway 

We consider the model proposed in [5] of a network responsible for the regulation of Apoptosis 
(cell death). This comprises the reactions listed below: 

C$ + C 3 
Cg + C3 
C* + IAP 
C£ + IAP 
C£ + BAR 





Due to the presence of many inflows as well as outflows (degradation of molecules) , the network 
does not exhibit conserved moieties; in the language of Petri Nets, there are no P-semiflows. 
Computation of T-semiflows, yields besides the obvious ones, R3 + R4, Rs + R9, R13 + Ru, 
i?i5 + Riq, Ru + R12, R17 + -Sis, entailed by the presence of reversible reactions, the following 
vectors: 



ct+ct 

CI -IAP -> 

CI -> 

C*- BAR - 

IAP 

BAR -» 

C 8 -> 

C 3 -» 

Ct - 



i?6 + Rll 
R2 + -^15 + R\9 
R\ + i?7 + i?i7 
i?l + i?3 + i?5 + -Ru + i?i7 
i?2 + i?8 + R\0 + -Rl3 + -Rl5 
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Notice that every reaction belongs to at least one T-semiflow, hence the network is consistent, 
and fulfills the necessary condition for persistence. Thus, it is interesting to look for sufficient 
conditions for persistence. Analysis of Input-Output stability of the network is made possible 
by the presence of outflows in every chemical species. Indeed, by letting 



V = [C 8 ] + [C$] + [C 3 ] + [C*] + [IAP] 

+ [BAR] + [C* - BAR] + [C* - IAP] 



(27) 



and taking its derivative with respect to time yields: 



V < (k 17 + k 15 + k 13 + k u )- k 16 [C 8 ] - fc 19 [C£] 
-k 18 [C 3 ] - k 7 [C$] - k 12 [IAP] - k u [BAR] 
-k 8 [C£ - BAR] - k 5 [C* - IAP] < K in - K out V 



(28) 



where the last inequality follows by letting: 



K, 



out 



= (k 17 + k 15 + k 13 + kn) 

= min{ fci 6 , feig , feis , k 7 , k 12 ,ku,k 8 ,k 5 }. 



(29) 



Hence, limsup t ^ +00 V(t) < Ki n /K out < +oo which proves input-output stability of the net- 
work (this is actually true even for time- varying inflows, in which case the sup norm of 
Ki n (t) is needed in constructing the estimate). To further analyze the network it is use- 
ful to investigate the presence of siphons. Indeed, there exists one minimal siphon: S = 
{C 3 ,C£,C 3 — IAP,C$ — BAR}. Due to the absence of P-semiflows this is indeed a critical 
siphon. Let a = [C 3 , C§ , C 3 — IAP, C§ — BAR]' be the concentration vector of species belonging 
to E and ( = [IAP, C 8 , C 3 , BAR]'. Writing down explicit equations for & yields: 



a 



-kjcri - fc 3 CTiCi + kiC, 3 a 2 + k 4 a 3 
-ki 9 a 2 - k 8 a 2 ( 4 + kga^ + k 2 ( 2 ai 
-(k 4 + &5)cr 3 + fc 3 <7iCi 

-(fcg + k W )(T4 + k 8 (J 2 (4 



(30) 



Taking the Jacobian of the vector-field for a = 0, yields: 



J 



-k 7 - k 3 (i 
k 2 ( 2 





h( 3 k A 

-fci 9 - k 8 (4 

— ki — k$ 

ha o 





kg 



-kg - fcio 



The species not included in E, on the other hand, evolve according to the following simple 
equations (once each variable in E is set to 0): 



Ci 


= k u 


~ &12& 


- k 3 &i(i - foviCi + ^40-3 




= k 15 


- k ie ( 2 


- k 2 a 2 Q 2 




= kyj 


- hs(3 


- ha 2 ( 3 


a 


= k 13 


- ku(4 


- k 8 a 2 (4 + kgG± 



(31) 
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which letting a = yields the following simple globally asymptotically stable embedded system: 



Ci = hi - fci2& 

C2 = k 15 ~ heb 

C3 = h7 - h 8 (3 

U = ki 3 - ku(4 



(32) 



whose equilibrium is located at [kn/ku, ki$/kiQ, kn/kis, ^13/^14]. Hence, we have a boundary 
equilibrium at [kn/ki 2 , fas/hie, kn/kis, hz/ku, 0, 0, 0, 0]. 

Applying Lemma IA.4I from the Appendix leads to the following. 

Lemma 10.1 Let Q > be arbitrary. Then J is Hurwitz if and only if its determinant is 
positive. □ 

Proof. 

Necessity is trivial. To prove sufficiency we calculate the 4 leading principal minors: 
M l = -kr-hd 

M 2 = (k 7 + k 3 (i) (ki 9 + k 8 Q) -^2(2(3 
M 3 = -(k 4 + k 5 )(k 7 (k 19 + kg( 4 )- ^2(2(3) 
-k3hCi(ki9 + k 8 Ci) 
= -(k 4 + h)M2 + k 4 k 3 ( 1 (k ig + kg( 4 ) 

M 4 = -(fcg + MMs-MM^ + A^ + A^sCi)/^ 

It follows that: 

det(J) = M 4 > => M 3 < ^ M 2 > 0, 

and Mi < is immediately clear. Thus, by Theorem IA.4I in the Appendix, J is Hurwitz 
whenever its determinant is positive. 



Consequently, if the determinant of J is positive, then the CRN is not persistent. On the other 
hand, if the determinant of J is negative, then it follows from Theorem [5] that the CRN is 
persistent (notice that J is irreducible). 



11 Conclusions 

New checkable criteria for persistence of chemical reaction networks have been proposed, which 
apply even when kinetic rates are time-dependent. These "time-dependent" rates may represent 
inflows and outflows, as well as the effect of external inputs. Finally, the case when critical 
siphons are present is also studied. As an illustration, a hypoxia and an apoptosis network are 
analyzed. 
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A Appendix 



In this Appendix we state and prove a few results used in the proof of Theorem 

Lemma A.l Assume that g : M n x W 71 —> M. m is a C 1 vector function in some open set 
containing and that it maps (x,y) to g(x,y). 

Suppose that there is some c^> such that: 

1. g(x,0) = for all x > 0. 

2. cdg/dy(x, 0) S> at some S >0. 

Then cg(x, y) > for all x in a neighborhood of x and all y 3> in a neighborhood of y = 0. □ 

Proof. For each (x,y), consider m(t) = g(x,ty) as a function of i in [0, 1]. By the Fundamental 
Theorem of Calculus, and using g(x, 0) = 0: 

g(x, y) = g(x, 0) + / (dm/dt)(t)dt = F(x,y)y, 
Jo 

where F(x,y) = f {dg /dy)(x,ty)dt. 

Pick a convex neighborhood Q of (x, 0), where cdg/dy(x, y) ^> for all (x, y) in Q. By convexity, 
also cdg/dy(x,ty) » for all t in [0, 1] and all (x,y) in Q; therefore cF(x,y) ^> for all (x,y) 
in Q, and hence 05(2;, y) > if (x, y) is in Q and y S> 0. I 

Lemma A. 2 Suppose that C is a compact, invariant set of a dynamical system generated by 
a system x = F(x) defined on some set X in R n , and that xo is a point such that: 

1. For every compact, forward invariant subset D of C, xq belongs to D. 

2. For every x in C different from xq, xq does not belong to a(x). 

Then, C = {x }. □ 

Proof. Suppose that there is some x in C which is different from xq. Consider D = a(x). Since 
C is compact and invariant, D is nonempty, compact, and invariant. Since D is, in particular, 
forward invariant, it follows from 1 that xq € D. A contradiction is obtained because 2 asserts 
that xq does not belong to D. I 



Corollary A. 3 Suppose that C is a compact, invariant set of a dynamical system generated 
by a system x = F{x) defined on some set X in M n , and that xo is globally asymptotically 
stable. Then, C = {x }. □ 
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Proof. Denote the flow on C generated by the system by (f>(t,x). We verify properties 1 and 
2 from Lemma lA.21 Property 1 is clear because of attractivity of xq. Property 2 follows from 
stability of xq. Indeed, pick a neighborhood U of xq which doesn't contain x, and a neighborhood 
V C U of xq such that trajectories starting in V cannot exit U in positive time. If xq belongs 
to a(x), then there is some t > such that <f>{—t,x) £ V. Then x = cfr(t,(j)(—t,x)) 6 U, a 
contradiction. I 

In the example of the regulation pathway for apoptosis, the following basic result from [6] is 
used. 

Lemma A. 4 A Metzler matrix is Hurwitz if and only if its leading principal minors alternate 
in sign, with the first one being negative. □ 
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